/**
  TADAQ Project -- online from TAAna Online by SUN Yazhou
  \class TAMath
  \brief Math class, to provide some general math methods.
  \author SUN Yazhou, aisa.rabbit@163.com
  \date Created: 2020/07/09
  \date Last modified: 2023/10/18 by SUN Yazhou
  \copyright 2020-2023 SUN Yazhou
  \copyright Copyright (c) 2021-2024 Anyang Normal U. && IMP-CAS with LGPLv3 LICENSE
  \copyright MOMD project, Anyang Normal University, IMP-CAS
*/

#include <iostream>
#include <cstring>
#include <algorithm>
#include <complex> // complex number operation
#include <utility>
#include <unistd.h> // for access
#include <cstdio>
#include <array>
#include <sys/stat.h>
#include <errno.h>
#include <ctype.h>

#include "TAMath.h"
#include "TAMsg.h"

using std::swap;
using std::max;
using std::cout;
using std::endl;
using std::array;

#define er TAMsg::Error

// returns the value ln[Gamma(xx)] for xx > 0.
double TAMath::gammaln(double xx){
  if(xx <= 0.) er("TAMath", "gammaln: x > 0. is mandatory.");
  // internal arithmetic will be done in double precision, a nicety that you can
  // omit if five-figure accuracy is good enough
  static const double c[6] = {76.18009172947146,-86.50532032941677,
    24.01409824083091,-1.231739572450155,
    0.1208650973866179e-2,-0.5395239384953e-5};
  double tmp = xx + 5.5, x = xx; // (z+gamma+1/2), gamma = 5, N = 6
  tmp -= (x+0.5)*log(tmp);
  double ser = 1.000000000190015; // the additive series
  for(int i = 0; i < 6; i++) ser += c[i]/++x;
  return -tmp+log(2.5066282746310005*ser/xx); // sqrt(2pi)
} // end member function gammaln

// returns ln(n!)
double TAMath::factln(int n){
  static double a[101]; // a static array is automatically initialized to zero

  if(n < 0) er("TAMath", "factln: n is negative.");
  if(n <= 1) return 0.;
  if(n <= 100) return a[n] ? a[n] : (a[n]=gammaln(n+1.));
  return gammaln(n+1.); // out of range of table
} // end member function factln

// \retval n!
double TAMath::Factorial(int n){
  static int ntop = 4;
  static const int NMAX = 101; // size of the cached table, max=171 for double
  static const int NMAX1 = NMAX - 1;
  static double a[NMAX] = {1., 1., 2., 6., 24.}; // fill in table only as required

  if(n < 0) er("TAMath", "Factorial: n: %d is minus", n);
  // larger value than size of the table is required. Actually, this big a value
  // is going to overflow on many computers, but no harm in trying
  if(n > NMAX1) return exp(gammaln(n+1.));
  while(ntop < n){
    ntop++;
    a[ntop] = a[ntop-1]*ntop;
  } // end while
  return floor(0.5+a[n]); // clear off the roundoff error
} // end member function Factorial

// returns Binomial coefficients
int TAMath::Binomial(int n, int m){
  if(m > n)
    er("TAMath",
      "Binomial: n: %d is larger than m: %d!", n, m);

  // the floor function cleans up roundoff error for smaller values of n and m
  return floor(0.5+exp(factln(n)-factln(m)-factln(n-m)));
} // end member function Binomial

double TAMath::Beta(double z, double w){
  if(z < 0. || w < 0.) er("TAMath", "Beta: z or z is minus");
  return exp(gammaln(z)+gammaln(w)-gammaln(z+w));
} // end member function Beta

// \retval n!!
int TAMath::BiFactorial(int n){
  return n <= 1 ? 1 : n * BiFactorial(n-2);
} // end member function BiFactorial

// the incomplete gamma function gammap, gammap(a,x)=P(a,x)=gamma(a,x)/Gamma(a)
double TAMath::gammap(double a, double x){
  if(x < 0. || a <= 0.) er("TAMath", "gammap: invalid arguements");

  if(x < (a+1.)) return gser(a, x); // series representation
  return 1. - gcf(a, x);
} // end member function gammap
// the incomplete gamma function gammaq, gammaq(a,x)=Q(a,x)=1-P(a,x)
double TAMath::gammaq(double a, double x){
  if(x < 0. || a <= 0.) er("TAMath", "gammaq: invalid arguements");

  if(x < (a+1.)) return 1. - gser(a, x); // continued fraction representation
  return gcf(a, x);
} // end member function gammaq

// returns the incomplete gamma function P(a,x) evaluated by its series
// representation. Optionally returns ln[Gamma(a)] as gln.
double TAMath::gser(double a, double x, double *gln){
  static const double ITMAX = 100; // maximum iteration times
  static const double EPS = 3.e-7; // fractional accuracy

  if(x < 0.) er("TAMath", "gser: x < 0., x: %f", x);
  double ap = a, sum, item = sum = 1./a;
  for(int i = 0; i < ITMAX; i++){
    sum += (item *= x/++ap);
    if(fabs(item) < fabs(sum)*EPS)
      return sum * exp(-x + a*log(x) - (gln ? *gln=gammaln(a) : gammaln(a)));
  } // end for over i
  er("TAMath", "gser: a too large, and IMAX too small.");
  return 0.; // never gets here
} // end member function gamser

// returns the incomplete gamma function Q(a,x) evaluated by its continued
// fraction representation. Optionally returns ln[Gamma(a)] as gln.
double TAMath::gcf(double a, double x, double *gln){
  static const double ITMAX = 100; // maximum iteration times
  static const double EPS = 3.e-7; // factional accuracy
  // number near the smallest representable floating-point number
  static const double FPMIN = 1.e-30; // modified Lentz's algorithm

  // set up for evaluating continued fraction by modified Lentz's method with b0=0.
  // b1 = x+1.-a, c1 = \infty (since c0 = 0), d1 = 1/b1, fr1=b0+a1/b1=a1/b1=1/b1=d1
  double b = x+1.-a, c = 1./FPMIN, d = 1./b, fr = d;
  for(int i = 1; i < ITMAX; i++){
    double aip = -i*(i-a); // a_{i+1}
    b += 2.; d = aip*d+b; // b2, d2 = 1/(a2*d1+b2)
    if(fabs(d) < FPMIN) d = FPMIN;
    c = b+aip/c; // c2 = b2 + a2/c1
    if(fabs(c) < FPMIN) c = FPMIN;
    d = 1./d;
    double delta = c*d;
    fr *= delta; // fr_2 = fr_1*c2*d2
    if(fabs(delta-1.) < EPS)
      return exp(-x+a*log(x)-(gln ? *gln = gammaln(a) : gammaln(a)))*fr; // put factors in front
  } // end for over i
  er("TAMath", "gcf: a too large, and IMAX too small.");
  return 0.; // never gets here
} // end member function gcf


double TAMath::BetaGamma(double beta){
  if(beta < 0. || beta >= 1.){
    // er("TAMath", "BetaGamma: input beta invalid: %f", beta);
    return -999.;
  } // end if
  return beta/sqrt(1.-beta*beta);
} // end member function BetaGamma

// linear interpolation
// kk stores the slope upon request
double TAMath::lii(double x0, double y0, double x1, double y1, double x, double *kk){
  // if(fabs(x0-x1) < fabs(y0-y1)*1.e-15) er("TAMath", "lii: the denominator too small");
  const double k = (y1-y0)/(x1-x0);
  if(kk) *kk = k;
  if(fabs(x - x0) < fabs(x - x1)) return k*(x-x0) + y0;
  else return k*(x-x1) + y1;
} // end member function lii

long TAMath::filesize(const string &f){
  if(access(f.data(), F_OK)){
    TAMsg::Warn("TAMath", "filesize: %s does not exist", f.data());
    return 0;
  } // end if
  static struct stat buf;
  if(-1 == stat(f.data(), &buf))
    er("TAMath", "filsize: %s - %s", f.data(), strerror(errno));
  return buf.st_size;
} // end member function filesize

double TAMath::filesizeInMB(const string &f){
  if(access(f.data(), F_OK)){
    TAMsg::Warn("TAMath", "filesizeInMB: %s does not exist", f.data());
    return 0;
  } // end if
  return filesize(f) / MB();
} // end member function filesize

bool TAMath::istrcmp(const string &s0, const string &s1){
  if(s0.size() != s1.size()) return false;
  for(int i = s0.size(); i--;)
    if(tolower(s0[i]) != tolower(s1[i])) return false;
  return true;
} // end member function istrcmp

void TAMath::uppercase(string &s){
  for(auto &c : s) c = toupper(c);
} // end member function toupper
